social_census_dashboard – Reddit Analysis of Cecil B. Moore Library and Carousel House

Code

import geopandas as gpd

#PASDA libraries shapefile
libraries = gpd.read_file("/Users/macbook/Desktop/Python/Week 14/PhiladelphiaLibraries201302/PhiladelphiaLibraries201302.shp")

#Swiming Pools 
pools = gpd.read_file("/Users/macbook/Desktop/Python/Week 14/PPR_Swimming_Pools/PPR_Swimming_Pools.shp")

#Playgrounds and Rec Centers
playgrounds = gpd.read_file("/Users/macbook/Desktop/Python/Week 14/PPR_Playgrounds/PPR_Playgrounds.shp")

#Parks
parks = gpd.read_file("/Users/macbook/Desktop/Python/Week 14/PPR_Program_Sites.geojson")

Code

for name, gdf in {
    "libraries": libraries,
    "pools": pools,
    "playgrounds": playgrounds,
    "parks": parks
}.items():
    print(name, gdf.crs)

libraries EPSG:2272
pools EPSG:3857
playgrounds EPSG:3857
parks EPSG:4326

Code

TARGET_CRS = "EPSG:2272"

layers = {
    "libraries": libraries,
    "pools": pools,
    "playgrounds": playgrounds,
    "parks": parks
}

for name, gdf in layers.items():
    if gdf.crs != TARGET_CRS:
        layers[name] = gdf.to_crs(TARGET_CRS)

libraries = layers["libraries"]
pools = layers["pools"]
playgrounds = layers["playgrounds"]
parks = layers["parks"]

Code

libraries["type"] = "Library"
pools["type"] = "Pool"
playgrounds["type"] = "Playground"
parks["type"] = "Park"

Code

import pandas as pd
import geopandas as gpd
import folium
import censusdata

Code

u18_vars = [
    "B01001_001E",  # total population
    "B01001_003E","B01001_004E","B01001_005E","B01001_006E",
    "B01001_027E","B01001_028E","B01001_029E","B01001_030E",
]

acs_u18 = censusdata.download(
    "acs5",
    2022,
    censusdata.censusgeo([("state","42"), ("county","101"), ("tract","*")]),
    u18_vars
).reset_index()

Code

edu_vars = [
    "B15003_001E",  # total 25+
    "B15003_022E","B15003_023E","B15003_024E","B15003_025E"
]

acs_edu = censusdata.download(
    "acs5",
    2022,
    censusdata.censusgeo([("state","42"), ("county","101"), ("tract","*")]),
    edu_vars
).reset_index()

Code

inc_vars = ["B19013_001E"]

acs_inc = censusdata.download(
    "acs5",
    2022,
    censusdata.censusgeo([("state","42"), ("county","101"), ("tract","*")]),
    inc_vars
).reset_index()

Code

def make_geoid(df):
    df["GEOID"] = df["index"].apply(
        lambda x: "".join([part[1] for part in x.geo])
    )
    return df

Code

acs_u18 = make_geoid(acs_u18)
acs_edu = make_geoid(acs_edu)
acs_inc = make_geoid(acs_inc)

Code

acs_u18["under_18"] = acs_u18[
    [
        "B01001_003E","B01001_004E","B01001_005E","B01001_006E",
        "B01001_027E","B01001_028E","B01001_029E","B01001_030E"
    ]
].sum(axis=1)

acs_u18["pct_under_18"] = (
    acs_u18["under_18"] / acs_u18["B01001_001E"] * 100
)

Code

acs_edu["bachelors_plus"] = acs_edu[
    ["B15003_022E","B15003_023E","B15003_024E","B15003_025E"]
].sum(axis=1)

acs_edu["pct_bachelors_plus"] = (
    acs_edu["bachelors_plus"] / acs_edu["B15003_001E"] * 100
)

Code

acs_inc = acs_inc.rename(columns={
    "B19013_001E": "median_income"
})

Code

acs_full = (
    acs_u18[["GEOID", "pct_under_18"]]
    .merge(
        acs_edu[["GEOID", "pct_bachelors_plus"]],
        on="GEOID",
        how="left"
    )
    .merge(
        acs_inc[["GEOID", "median_income"]],
        on="GEOID",
        how="left"
    )
)

Code

import geopandas as gpd

# Read tracts geometry and merge the consolidated ACS table (acs_full)
tracts_geom = gpd.read_file(
    "https://www2.census.gov/geo/tiger/TIGER2022/TRACT/tl_2022_42_tract.zip"
)

tracts_geom = tracts_geom[tracts_geom["COUNTYFP"] == "101"]

tracts = tracts_geom.merge(
    acs_full,
    on="GEOID",
    how="left",
)

Code

# --- social infrastructure attribute renames and combine ---
libraries = libraries.rename(columns={"ASSET_NAME": "name"})
pools = pools.rename(columns={"amenity_na": "name"})
playgrounds = playgrounds.rename(columns={"park_name": "name"})
parks = parks.rename(columns={"park_name": "name"})

Code

social_infra = gpd.GeoDataFrame(
    pd.concat([
        libraries[["name", "type", "geometry"]],
        pools[["name", "type", "geometry"]],
        playgrounds[["name", "type", "geometry"]],
        parks[["name", "type", "geometry"]],
    ], ignore_index=True),
    crs=libraries.crs
)

Code

social_infra_web = social_infra.to_crs("EPSG:4326")

Code

# Project tracts and points to a projected CRS
tracts_proj = tracts.to_crs("EPSG:2272")
social_infra_proj = social_infra.to_crs("EPSG:2272")

Code

infra_join = gpd.sjoin(
    social_infra_proj,
    tracts_proj[["GEOID", "geometry"]],
    how="left",
    predicate="within"
)

Code

infra_counts = (
    infra_join
    .groupby("GEOID")
    .size()
    .reset_index(name="facility_count")
)

Code

infra_counts.head()
infra_counts["facility_count"].describe()

count    225.000000
mean       3.382222
std        2.032250
min        1.000000
25%        2.000000
50%        3.000000
75%        4.000000
max       13.000000
Name: facility_count, dtype: float64

Code

type_cols = {
    "Library": "library_count",
    "Playground": "rec_center_count",
    "Rec Center": "rec_center_count",
    "Park": "park_count",
    "Pool": "pool_count"
}
infra_by_type = (
    infra_join
    .groupby(["GEOID", "type"])
    .size()
    .unstack(fill_value=0)
    .reset_index()
    .rename(columns=type_cols)
)
for col in type_cols.values():
    if col not in infra_by_type.columns:
        infra_by_type[col] = 0

Code

tracts_access = (
    tracts_proj
    .merge(infra_counts, on="GEOID", how="left")
    .merge(infra_by_type, on="GEOID", how="left")
)

# Replace NA with 0
count_cols = ["facility_count", "library_count", "rec_center_count", "park_count", "pool_count"]
tracts_access[count_cols] = tracts_access[count_cols].fillna(0)

Code

tracts_access_web = tracts_access.to_crs("EPSG:4326")

Code

count_field_options = [
    ("Social Infrastructure Count:", ["facility_count"]),
    ("Libraries:", ["library_count", "Library"]),
    ("Rec Centers:", ["rec_center_count", "Rec Center", "Playground"]),
    ("Parks:", ["park_count", "Park"]),
    ("Pools:", ["pool_count", "Pool"]),
]
tooltip_fields = []
tooltip_aliases = []
for alias, options in count_field_options:
    for col in options:
        if col in tracts_access_web.columns:
            tooltip_fields.append(col)
            tooltip_aliases.append(alias)
            break
if not tooltip_fields:
    tooltip_fields = ["facility_count"]
    tooltip_aliases = ["Social Infrastructure Count:"]

Code

from folium.plugins import DualMap

dual_map = DualMap(
    location=[39.9526, -75.1652],
    zoom_start=11,
    tiles="CartoDB positron"
)

Code

import numpy as np

def style_facilities(feature):
    v = feature["properties"].get("facility_count")
    if v is None or v == 0 or (isinstance(v, float) and np.isnan(v)):
        return {
            "fillColor": "#f0f0f0",
            "color": "#cccccc",
            "weight": 0.3,
            "fillOpacity": 0.45
        }
    return {
        "fillColor": cmap_facilities(v),
        "color": "#ffffff",
        "weight": 0.3,
        "fillOpacity": 0.9
    }

Code

tracts_access_web["facility_count"].describe()

count    408.000000
mean       1.865196
std        2.260431
min        0.000000
25%        0.000000
50%        1.000000
75%        3.000000
max       13.000000
Name: facility_count, dtype: float64

Code

import branca.colormap as cm


bins = tracts_access_web["facility_count"].quantile(
    [0, 0.2, 0.4, 0.6, 0.8, 1]
).values

cmap_facilities = cm.StepColormap(
    colors=cm.linear.Blues_06.colors,
    index=bins,
    vmin=bins[0],
    vmax=bins[-1]
)

Code

import censusdata
import geopandas as gpd
import pandas as pd
import folium

Code

acs_vars = {
    "total_pop": "B01001_001E",
    "male_under_18": [
        "B01001_003E","B01001_004E","B01001_005E","B01001_006E"
    ],
    "female_under_18": [
        "B01001_027E","B01001_028E","B01001_029E","B01001_030E"
    ]
}

tracts_acs = censusdata.download(
    "acs5",
    2022,
    censusdata.censusgeo([("state","42"), ("county","101"), ("tract","*")]),
    [acs_vars["total_pop"]] +
    acs_vars["male_under_18"] +
    acs_vars["female_under_18"]
)

Code

tracts_acs["under_18"] = tracts_acs[
    [
        "B01001_003E",
        "B01001_004E",
        "B01001_005E",
        "B01001_006E",
        "B01001_027E",
        "B01001_028E",
        "B01001_029E",
        "B01001_030E",
    ]
].sum(axis=1)

Code

edu_vars = [
    "B15003_001E",  # total 25+
    "B15003_022E",  # bachelor's
    "B15003_023E",  # master's
    "B15003_024E",  # professional
    "B15003_025E"   # doctorate
]

Code

income_vars = ["B19013_001E"]

Code

tracts_acs_edu = censusdata.download(
    "acs5",
    2022,
    censusdata.censusgeo([("state","42"), ("county","101"), ("tract","*")]),
    edu_vars
)

tracts_acs_inc = censusdata.download(
    "acs5",
    2022,
    censusdata.censusgeo([("state","42"), ("county","101"), ("tract","*")]),
    income_vars
)

Code

tracts_acs_edu = tracts_acs_edu.reset_index()
tracts_acs_inc = tracts_acs_inc.reset_index()

Code

tracts_acs_edu = tracts_acs_edu.rename(columns={
    "B15003_001E": "total_25plus",
    "B15003_022E": "bachelors",
    "B15003_023E": "masters",
    "B15003_024E": "professional",
    "B15003_025E": "doctorate"
})

Code

# If total population is still B01001_001E, rename it first
tracts_acs = tracts_acs.rename(columns={
    "B01001_001E": "total_pop"
})

# Recreate under_18 if needed
if "under_18" not in tracts_acs.columns:
    tracts_acs["under_18"] = tracts_acs[
        [
            "B01001_003E","B01001_004E","B01001_005E","B01001_006E",
            "B01001_027E","B01001_028E","B01001_029E","B01001_030E"
        ]
    ].sum(axis=1)

# NOW create the percentage
tracts_acs["pct_under_18"] = (
    tracts_acs["under_18"] / tracts_acs["total_pop"] * 100
)

Code


if "bachelors" in tracts_acs_edu.columns and "pct_bachelors" not in tracts_acs_edu.columns:
    tracts_acs_edu["pct_bachelors"] = (
        tracts_acs_edu["bachelors"]
        / tracts_acs_edu["total_25plus"] * 100
    )

Code

tracts_acs_inc = tracts_acs_inc.rename(columns={
    "B19013_001E": "median_income"
})

Code

tracts_acs_edu["GEOID"] = tracts_acs_edu["index"].apply(
    lambda x: "".join([part[1] for part in x.geo])
)

Code

tracts_acs_inc["GEOID"] = tracts_acs_inc["index"].apply(
    lambda x: "".join([part[1] for part in x.geo])
)

Code

tracts_acs["GEOID"] = tracts_acs.index.map(
    lambda x: (
        x.geo[0][1] +  # state
        x.geo[1][1] +  # county
        x.geo[2][1]    # tract
    )
)

Code

tracts_acs[["total_pop", "under_18", "pct_under_18"]].head()

	total_pop	under_18	pct_under_18
Census Tract 1.01; Philadelphia County; Pennsylvania: Summary level: 140, state:42> county:101> tract:000101	1947	40	2.054443
Census Tract 1.02; Philadelphia County; Pennsylvania: Summary level: 140, state:42> county:101> tract:000102	2897	162	5.591992
Census Tract 2; Philadelphia County; Pennsylvania: Summary level: 140, state:42> county:101> tract:000200	3486	414	11.876076
Census Tract 3; Philadelphia County; Pennsylvania: Summary level: 140, state:42> county:101> tract:000300	3914	280	7.153807
Census Tract 4.01; Philadelphia County; Pennsylvania: Summary level: 140, state:42> county:101> tract:000401	2675	127	4.747664

Code

acs_full = (
    tracts_acs[["GEOID", "pct_under_18"]]
    .merge(
        tracts_acs_edu[["GEOID", "pct_bachelors"]],
        on="GEOID",
        how="left"
    )
    .merge(
        tracts_acs_inc[["GEOID", "median_income"]],
        on="GEOID",
        how="left"
    )
)

Code

tracts = tracts_geom.merge(
    acs_full,
    on="GEOID",
    how="left"
)

tracts_web = tracts.to_crs("EPSG:4326")

Code

import folium

m_acs = folium.Map(
    location=[39.9526, -75.1652],
    zoom_start=11,
    tiles="CartoDB positron"
)

Code

import numpy as np

bins_income = tracts_web["median_income"].quantile(
    [0, 0.2, 0.4, 0.6, 0.8, 1]
).values

Code

import branca.colormap as cm

cmap_income = cm.StepColormap(
    colors=cm.linear.YlOrRd_06.colors,
    vmin=bins_income[0],
    vmax=bins_income[-1],
    index=bins_income
)

Code

bins_bach = tracts_access_web["pct_bachelors_plus"].quantile(
    [0, .2, .4, .6, .8, 1]
).values

Code

import branca.colormap as cm
import numpy as np

cmap_bach = cm.StepColormap(
    colors=cm.linear.Purples_06.colors,
    index=bins_bach,
    vmin=bins_bach[0],
    vmax=bins_bach[-1]
)
cmap_bach.caption = "% Bachelor's Degree or Higher"

Code

def style_bachelors(feature):
    v = feature["properties"].get("pct_bachelors_plus")
    if v is None or np.isnan(v):
        return {
            "fillColor": "#A9A9A9",
            "color": "white",
            "weight": 0.4,
            "fillOpacity": 0.85
        }
    return {
        "fillColor": cmap_bach(v),
        "color": "white",
        "weight": 0.4,
        "fillOpacity": 0.85
    }

Code

folium.GeoJson(
    tracts_web,
    name="% Bachelor’s Degree or Higher",
    style_function=style_bachelors
).add_to(m_acs)

<folium.features.GeoJson at 0x14ee4e5d0>

Code

cmap_income = cm.linear.YlOrRd_09.scale(
    tracts_web["median_income"].min(),
    tracts_web["median_income"].max()
)

Code

import numpy as np

bins_income = tracts_web["median_income"].quantile(
    [0, 0.2, 0.4, 0.6, 0.8, 1]
).values

Code

import branca.colormap as cm

cmap_income = cm.StepColormap(
    colors=cm.linear.YlOrRd_06.colors,
    vmin=bins_income[0],
    vmax=bins_income[-1],
    index=bins_income
)

Code

def build_acs_layer(name, field, palette):
    values = tracts_access_web[field]
    bins = values.quantile([0, 0.2, 0.4, 0.6, 0.8, 1]).values
    cmap = cm.StepColormap(
        colors=palette,
        index=bins,
        vmin=bins[0],
        vmax=bins[-1]
    )
    cmap.caption = name
    def style_fn(feature, field=field, cmap=cmap):
        val = feature["properties"].get(field)
        if val is None or (isinstance(val, float) and np.isnan(val)):
            return {"fillColor": "#d9d9d9", "color": "white", "weight": 0.3, "fillOpacity": 0.4}
        return {"fillColor": cmap(val), "color": "white", "weight": 0.3, "fillOpacity": 0.85}
    return {"name": name, "field": field, "style": style_fn, "cmap": cmap}

acs_layers = [
    build_acs_layer("% Bachelor's Degree or Higher", "pct_bachelors_plus", cm.linear.Purples_06.colors),
    build_acs_layer("% Under Age 18", "pct_under_18", cm.linear.Greens_06.colors),
    build_acs_layer("Median Household Income", "median_income", cm.linear.YlOrRd_06.colors),
]

FINAL MAP

Code

import folium
from folium.plugins import DualMap
from folium.features import GeoJsonTooltip
import branca.colormap as cm
import numpy as np

Code

#Dual Map Setup

from folium.plugins import DualMap
import folium
from folium.features import GeoJsonTooltip

dual_map_edu = DualMap(
    location=[39.9526, -75.1652],
    zoom_start=11,
    tiles="CartoDB positron"
)

Code

folium.GeoJson(
    tracts_access_web,
    style_function=style_facilities,
    tooltip=GeoJsonTooltip(
        fields=tooltip_fields,
        aliases=tooltip_aliases,
        localize=True,
        sticky=True
    )
).add_to(dual_map_edu.m1)

# Legend (optional but fine here)
cmap_facilities.add_to(dual_map_edu.m1)

Code

for layer in acs_layers:
    fg = folium.FeatureGroup(name=layer["name"], show=(layer["field"] == "pct_bachelors_plus"))
    folium.GeoJson(
        tracts_access_web,
        style_function=layer["style"],
        tooltip=GeoJsonTooltip(
            fields=[layer["field"]],
            aliases=[f"{layer['name']}:"],
            localize=True,
            sticky=True
        )
    ).add_to(fg)
    fg.add_to(dual_map_edu.m2)
    layer["cmap"].add_to(dual_map_edu.m2)

folium.LayerControl(collapsed=False).add_to(dual_map_edu.m2)

Code

tracts_access_web.head()

	STATEFP	COUNTYFP	TRACTCE	GEOID	NAME	NAMELSAD	MTFCC	FUNCSTAT	ALAND	AWATER	...	INTPTLON	geometry	pct_under_18	pct_bachelors_plus	median_income	facility_count	Park	Playground	Pool
0	42	101	026301	42101026301	263.01	Census Tract 263.01	G5020	S	406340	0	...	-075.1637161	POLYGON ((-75.16899 40.07147, -75.16864 40.071...	22.435737	28.359827	57992	0.0	0.0	0.0	0.0
1	42	101	029200	42101029200	292	Census Tract 292	G5020	S	1745920	28444	...	-075.1017831	POLYGON ((-75.11288 40.02746, -75.1128 40.0274...	23.318995	19.735431	82784	0.0	0.0	0.0	0.0
2	42	101	024400	42101024400	244	Census Tract 244	G5020	S	421189	0	...	-075.1638925	POLYGON ((-75.16912 40.02387, -75.16843 40.024...	33.115265	40.127077	56912	3.0	1.0	2.0	0.0
3	42	101	033200	42101033200	332	Census Tract 332	G5020	S	842716	0	...	-075.0449684	POLYGON ((-75.05464 40.04465, -75.05396 40.045...	29.119272	26.290400	80409	2.0	1.0	0.0	1.0
4	42	101	980200	42101980200	9802	Census Tract 9802	G5020	S	4916641	250846	...	-075.0443487	POLYGON ((-75.07448 40.08685, -75.07443 40.087...	6.666667	24.867725	-666666666	5.0	1.0	4.0	0.0

5 rows × 21 columns

Code

dual_map_edu

Make this Notebook Trusted to load map: File -> Trust Notebook

Code

folium.GeoJson(
    tracts_access_web,
    style_function=style_facilities,
    tooltip=GeoJsonTooltip(
        fields=tooltip_fields,
        aliases=tooltip_aliases,
        localize=True,
        sticky=True
    )
).add_to(dual_map.m1)

<folium.features.GeoJson at 0x16ba98150>

Code

for layer in acs_layers:
    fg = folium.FeatureGroup(name=layer["name"], show=(layer["field"] == "pct_bachelors_plus"))
    folium.GeoJson(
        tracts_access_web,
        style_function=layer["style"],
        tooltip=GeoJsonTooltip(
            fields=[layer["field"]],
            aliases=[f"{layer['name']}:"],
            localize=True,
            sticky=True
        )
    ).add_to(fg)
    fg.add_to(dual_map.m2)
    layer["cmap"].add_to(dual_map.m2)

folium.LayerControl(collapsed=False).add_to(dual_map.m2)

<folium.features.GeoJson at 0x16b7a7d10>

Code

# Add aligned legends for both map panes
cmap_facilities.add_to(dual_map.m1)

Code

dual_map

Make this Notebook Trusted to load map: File -> Trust Notebook